
im = dblarr(260,40,400)
mkhdr,hd,im
sz = size(im,/dim)



;; Create FITS astrometry structure
make_astr,astr,CRPIX=[sz[0]+0.5,sz[1]/2.+0.5],CRVAL=[0.,0.],$
          DELT=[-0.25d,0.05d],CTYPE=['GLON-CAR','GLAT-CAR']
putast,hd,astr

;; Add FITS header keywords and check NAXIS* count
check_fits,im,hd,/update
sxaddpar,hd,'BUNIT','f_fore',' Foreground fraction of IRAC Band 4 emission',$
         AFTER='EQUINOX'


;; Add information about the 3rd dimension...
sxaddpar,hd,'CDELT3',50.d,' Parsecs / Pixel',AFTER='CD2_2'
sxaddpar,hd,'CRPIX3',1.d,' Reference Pixel in Z',AFTER='CRPIX2'
sxaddpar,hd,'CRVAL3',50.d,' Heliocentric distance of reference pixel',$
         AFTER='CRVAL2'
sxaddpar,hd,'CTYPE3','DIST-HEL',' Coordinate Type',AFTER='CTYPE2'
sxaddpar,hd,'CUNIT1','DEGREE',' Axis 1 Units',AFTER='CTYPE1'
sxaddpar,hd,'CUNIT2','DEGREE',' Axis 2 Units',AFTER='CTYPE2'
sxaddpar,hd,'CUNIT3','PARSEC',' Axis 3 Units',AFTER='CTYPE3'

;; Write image + header to file
print,hd









cubed = dindgen(400)/400.*20000. + 50.
print,m4_stat(cubed)


for xi=0,259 do begin
   print,'Working column',xi
   for yi=0,39 do begin

      xy2ad,xi,yi,astr,l,b
      im[xi,yi,*] = ffore_lbd(l,b,cubed)
      
   endfor
endfor

print,m4_stat(im)

writefits,'sample_cube.fits',im,hd

END
